{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load all the required packages\n",
    "using ODEInterface\n",
    "using ForwardDiff\n",
    "using Gadfly\n",
    "using Colors\n",
    "@ODEInterface.import_huge\n",
    "loadODESolvers();\n",
    "\n",
    "# Define the right-hand function for Automatic Differentiation\n",
    "function roberAD(x)\n",
    "    return [-0.04*x[1]+1e4*x[2]*x[3],\n",
    "        0.04*x[1]-1e4*x[2]*x[3]-3e7*(x[2])^2,\n",
    "        3*10^7*(x[2])^2]\n",
    "end\n",
    "\n",
    "# Define the system for the solver\n",
    "function rober(t,x,dx)\n",
    "    dx[1] = -0.04*x[1]+1e4*x[2]*x[3];\n",
    "    dx[2] = 0.04*x[1]-1e4*x[2]*x[3]-3e7*(x[2])^2;\n",
    "    dx[3] = 3e7*(x[2])^2;\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# Automatic Differentiation for a more general problem\n",
    "function getJacobian(t,x,J)\n",
    "    J[:,:] = ForwardDiff.jacobian(roberAD,x);\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# Flag to check whether plot is to be generated and saved or not\n",
    "# Also checks if all solvers are successful\n",
    "printFlag = true;\n",
    "\n",
    "# Initial conditions\n",
    "t0 = 0.0; T = 10.^[0.0:11.0;]; x0=[1.0,0.0,0.0];\n",
    "\n",
    "# Get \"reference solution\"\n",
    "# TolMin < 1e-14 gives error \"TOLERANCES ARE TOO SMALL\" \n",
    "Tol = 1e-14;\n",
    "opt = OptionsODE(OPT_RHS_CALLMODE => RHS_CALL_INSITU,\n",
    "    OPT_RTOL => Tol, OPT_ATOL=>Tol*1e-6,OPT_EPS => 1.11e-16,\n",
    "    OPT_JACOBIMATRIX => getJacobian);\n",
    "\n",
    "(t,x,retcode,stats) = odecall(seulex,rober,[t0, T[end]], x0, opt);\n",
    "t = [t0;t];\n",
    "x = [x0';x];\n",
    "\n",
    "plotColorsHex = [\"#4D4D4D\",\"#5DA5DA\",\"#FAA43A\"];\n",
    "plotColors = [parse(Colorant,c) for c in plotColorsHex];\n",
    "\n",
    "p = plot(layer(x=t,y=x[:,1],Geom.path,Theme(line_width=2pt,default_color=plotColors[1])),\n",
    "layer(x=t,y=1e4*x[:,2],Geom.path,Theme(line_width=2pt,default_color=plotColors[2])),\n",
    "layer(x=t,y=x[:,3],Geom.path,Theme(line_width=2pt,default_color=plotColors[3])),Scale.x_log10,\n",
    "Guide.xlabel(\"time (s)\"), Guide.ylabel(\"Species concentration\"),\n",
    "Guide.manual_color_key(\"Legend\",[\"Specie 1\", \"(Specie 2)*1e4\", \"Specie 3\" ],plotColorsHex[1:3]),\n",
    "Theme(panel_stroke=colorant\"black\",highlight_width=0.1pt,key_position=:top,\n",
    "key_max_columns = 1,major_label_font_size=14pt,minor_label_font_size=10pt,\n",
    "key_title_font_size=12pt,key_label_font_size=10pt),\n",
    "Guide.xticks(ticks = [11:-1:-8;]),Guide.yticks(ticks=[1.0:-0.1:0.0;]));\n",
    "\n",
    "draw(PNG(\"../../ImagesAndPDFs/Plots/RoberPlot.png\",20cm,15cm),p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
